da=read.table("w-petroprice.txt",header = T)
da1=read.table("w-gasoline.txt")
pus=log(da$US)
pgs=log(da1[,1])
tdx=c(1:717)/52+1997
par(mfcol=c(2,1))
plot(tdx,pgs,xlab='year',ylab='ln(price)',type='l')
title(main='(a) Gasline')
plot(tdx,pus,xlab='year',ylab='ln(price)',type='l')
title(main='(b) Crude oil')
dpgs=diff(pgs)
acf(dpgs,lag=20)
pacf(dpgs,lag=20)
m1=ar(diff(pgs),method = 'mle')
t.test(dpgs)
m1=arima(dpgs,order=c(5,0,0),include.mean = F)
m1=arima(dpgs,order=c(5,0,0),include.mean = F,fixed = c(NA,NA,NA,0,NA))
tsdiag(m1,gof=20)
dpus=diff(pus)
m3=lm(dpgs~1+dpus)
acf(m3$residuals,lag=20)
pacf(m3$residuals,lag=20)
m4=ar(m3$residuals,method='mle')
m4=arima(dpgs,order=c(6,0,0),include.mean = F,xreg=dpus)
m4=arima(dpgs,order=c(5,0,0),include.mean = F,xreg=dpus)
m4=arima(dpgs,order=c(5,0,0),include.mean = F,xreg=dpus,fixed = c(NA,NA,NA,0,NA,NA))
tsdiag(m4,gof=20)
c1=c(NA,NA,NA,0,NA)
source('backtest.R')
pm1=backtest(m1,dpgs,316,1,fixed = c1,inc.mean = F)
c4=c(NA,NA,NA,0,NA,NA)
pm4=backtest(m4,dpgs,316,1,xre=dpus,inc.mean = F,fixed = c4)
tdx=tdx[2:717]
pm4fit=dpgs[317:716]-pm4$error
pm1fit=dpgs[317:716]-pm1$error
plot(tdx[317:716],dpgs[317:716],xlab='year',ylab='growth',type='l')
points(tdx[317:716],pm1fit,pch='*')
plot(tdx[317:716],dpgs[317:716],xlab='year',ylab='growth',type='l')
points(tdx[317:716],pm4fit,pch='*')
m6=lm(dpgs[2:716]~-1+dpus[1:715])
acf(m6$residuals,lag=20)
pacf(m6$residuals,lag=20)
m7=ar(m6$residuals,method = 'mle')
m7=arima(dpgs[2:716],order = c(9,0,0),xreg=dpus[1:715],include.mean = F)
m7=arima(dpgs[2:716],order=c(9,0,0),xreg=dpus[1:715],include.mean = F,fixed=c(NA,NA,NA,0,NA,0,0,0,NA,NA))
tsdiag(m7,gof=20)
c7=c(NA,NA,NA,0,NA,0,0,0,NA,NA)
pm7=backtest(m7,dpgs[2:716],315,1,xre=dpus[1:715],inc.mean = F,fixed = c7)

Gt=scan(file='m-GLBTs.txt')
Gtemp=ts(Gt,frequency = 12,start = c(1880,1))
plot(Gtemp,xlab='year',ylab='temperature',type='l')
acf(diff(Gt),lag=36)
pacf(diff(Gt),lag=36)
m1=arima(Gt,order=c(1,1,2))
acf(m1$residuals,lag=36)
m1=arima(Gt,order = c(1,1,2),seasonal = list(order=c(0,0,1),period=24))
tsdiag(m1,gof=36)
time=c(1:1568)
m2=lm(Gt~time)
par(mfcol=c(2,1))
acf(m2$residuals,lag=36)
m2=arima(Gt,order = c(2,0,1),xreg=time)
tsdiag(m2,gof=36)
m2=arima(Gt,order=c(2,0,1),seasonal = list(order=c(0,0,1),period=24),xreg=time)
tsdiag(m2,gof=36)
source('backtest.R')
pm1=backtest(m1,Gt,1368,1)
time=as.matrix(time)
pm2=backtest(m2,Gt,1368,1,xre=time)
time=c(1:1568)
time1=c(rep(0,1212),time[1213:1568])
mm1=lm(Gt~time+time1)
x1=cbind(time,time1)
mm1=arima(Gt,order=c(2,0,1),seasonal = list(order=c(0,0,1),period=24),xreg=x1)
tsdiag(mm1,gof=36)
Box.test(mm1$residuals,lag=8,type='Ljung')

da=read.table('m-unrate.txt',header=T)
unemp=da$rate
unrate=ts(unemp,frequency = 12,start = c(1948,1))
plot(unrate,xlab='year',ylab='unrate',type='l')
par(mfcol=c(2,2))
acf(unemp,lag=36)
pacf(unemp,lag=36)
acf(diff(unemp),lag=36)
pacf(diff(unemp),lag=36)
m1=arima(unemp,order=c(1,1,5),seasonal = list(order=c(1,0,1),period=12))
c1=c(NA,NA,NA,0,0,NA,NA,NA)
m1=arima(unemp,order=c(1,1,5),seasonal = list(order=c(1,0,1),period=12),fixed=c1)
tsdiag(m1,gof=36)
Box.test(m1$residuals,lag=24,type='Ljung')
Box.test(m1$residuals,lag=36,type='Ljung')
mm=arima(unemp,order=c(0,1,0),seasonal = list(order=c(1,0,1),period=12))
par(mfcol=c(2,1))
acf(mm$residuals,lag=24)
pacf(mm$residuals,lag=24)
mm1=arima(unemp,order = c(5,1,0),seasonal = list(order=c(1,0,1),period=12))
cc1=c(0,NA,NA,NA,NA,NA,NA)
mm1=arima(unemp,order=c(5,1,0),seasonal = list(order=c(1,0,1),period=12),fixed = cc1)
tsdiag(mm1,gof=36)
source('backtest.R')
pm1=backtest(m1,unemp,700,1,fixed = c1,inc.mean = F)
pmm1=backtest(mm1,unemp,700,1,fixed = cc1,inc.mean = F)

da=read.table('m-unrateic.txt',header = T)
unrate=da$rate
x=da[,5:9]/1000
nm1=lm(unrate~icm1,data=x)
